#############################
#############################
## Replication file #########
## BJPS submission ##########
## Nov 15 2023  #############
#############################
#############################

rm(list = ls())

## Analysis conducted in base R version 4.3.0

## set wd

## load packages 

library(ggplot2)
library(directlabels)
library(RColorBrewer)
library(colorRamps)
getPalette = colorRampPalette(brewer.pal(9, "Spectral"))
library(lme4)
library(lmerTest)
library(interplot)
library(dplyr)
library(plyr)
library(reshape2)
library(broom)
library(purrr)
library(stargazer)
library(xtable)

### load data 

load("prolific_data_bjps_replication.RData")
head(dat)

#################################
### Codebook -- Variable names:
### woman = Woman 
### highered = Higher education
### white = White 
### Age = Age 
### full_time = Works full-time 
### age_youngest_1 = Age of youngest child
### no_child2 = Number of children living at home 
### ml_share = Mental load reported (personal) share
### ml_open_count = Number of characters, mental load response 
### ml_satisfied = Satisfied with mental load 
### ml_fair_yes = Fairness perception of mental load 
### sneg_emotions = Negative emotions about mental load
### treated = Treated (mental load priming treatment)
### sinterest = Political interest 
### spol_ambition_scale = Political ambition
### spol_vote_scale = Vote intention
### spol_public  = Public participation
### spol_private = Private participation
### work_scale1_2 = Workplace advancement
### work_hrs = (preferred) work hours 
### work_hrs_partner = (preferred) partner work hours

###############################
### Table 1, summary stats
###############################

vars<-c("woman", "Age", "white", "highered", "full_time", "age_youngest_1", "no_child2", "ml_share")
sumdat<-dat[,vars]
stargazer(sumdat)

###############################
### Figure 1
###############################

summary(dat$ml_share)

dat$Gender<-NA
dat$Gender[dat$woman==1]<-"Woman"
dat$Gender[dat$woman==0]<-"Man"

vars<-c("Gender", "ml_share")
temp<-na.omit(dat[,vars])
nrow(temp)

mu <- ddply(temp, "Gender", summarise, grp.mean=mean(ml_share))
head(mu)

mu$grp.mean[mu$Gender=="Woman"]<-mean(dat$ml_share[dat$woman==1], na.rm=TRUE)

p<- ggplot(temp, aes(x=ml_share, colour=Gender, fill=Gender)) + 
  geom_density(alpha=0.2) +
  geom_vline(data=mu, aes(xintercept=grp.mean, color=Gender),
             linetype="dashed") + theme_classic()+
labs(title="",x="Share of cognitive household labor (self-reported estimate of total done by me)", y = "Density") 

p 

ggsave("Fig1.tiff", plot = p, width = 6, height = 4, dpi = 600, compression = "lzw")

p_grey <- ggplot(temp, aes(x=ml_share, colour=Gender, fill=Gender)) + 
  geom_density(alpha=0.2) +
  geom_vline(data=mu, aes(xintercept=grp.mean, color=Gender),
             linetype="dashed") + 
  theme_classic() +
  labs(title="",x="Share of cognitive household labor (self-reported estimate of total done by me)", y = "Density") +
  scale_color_grey(start = 0.3, end = 0.7) +   # Convert color to greyscale
  scale_fill_grey(start = 0.3, end = 0.7)      # Convert fill to greyscale

ggsave("Fig1_greyscale.tiff", plot = p_grey, width = 6, height = 4, dpi = 600, compression = "lzw")


t9 <- t.test(ml_share~ Gender, data=temp)
t9

median(dat$ml_share[dat$woman==0], na.rm=TRUE)
median(dat$ml_share[dat$woman==1], na.rm=TRUE)

######################################
##### Table 2
##### Mental load questions by gender
#######################################

t1a <- t.test(ml_share~ woman, data=dat)
t2a <- t.test(ml_satisfied~ woman, data=dat)
t3a <- t.test(ml_fair_yes~ woman, data=dat)
t4a <- t.test(sneg_emotions ~ woman, data=dat)
t5a <- t.test(ml_open_count ~ woman, data=dat)

tab <- map_df(list(t1a, t5a, t2a, t3a, t4a), tidy)

gender_div<-tab[c("estimate1", "estimate2", "estimate", "p.value")]
xtable(gender_div, type = "latex")

### Does number of characters correlate with ml_share? 

cor(dat$ml_share, dat$ml_open_count,  method = "pearson", use = "complete.obs")
cor.test(dat$ml_share, dat$ml_open_count,  method = "pearson", use = "complete.obs") ## it is sig but fairly weak

### among women 
cor.test(dat$ml_share[dat$woman==1], dat$ml_open_count[dat$woman==1],  method = "pearson", use = "complete.obs") ## it is sig but fairly weak

### among men 
cor.test(dat$ml_share[dat$woman==0], dat$ml_open_count[dat$woman==0],  method = "pearson", use = "complete.obs") ## it is sig but fairly weak

#########################################
### Figure 3 (note: Figure 2 code is at end)
#########################################

mod1<-lm(sinterest~treated, data=dat)
summary(mod1)
mod2<-lm(spol_ambition_scale~treated, data=dat)
summary(mod2)
mod3<-lm(spol_vote_scale~treated, data=dat)
summary(mod3)
mod4<-lm(spol_public~treated, data=dat)
summary(mod4)
mod5<-lm(spol_private~treated, data=dat)
summary(mod5)
mod6<-lm(work_scale1_2~treated, data=dat)
summary(mod6)

results <- list(mod1, mod2, mod3, mod4, mod5, mod6)

fx <- sapply(results, function(x){unlist(summary(x)$coefficients["treated", 1])})
l95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 1])})
u95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 2])})
ys <- c(6,5,4,3,2,1)

tiff("Fig3.tiff", height = 6, width = 10, units = "in", res = 600, compression = "lzw")

par(mar=c(5.1,10,2.1,2.1))

plot(fx, ys, xlim=c(-.1,.05), ylab='', yaxt='n', type='n', xaxt='n', main="", xlab="")
abline(v=0)

for (i in 1:6){
  if(i %% 2 != 0){
    segments(l95[i], ys[i], u95[i], ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  } else {
    segments(l95[i], ys[i], u95[i], ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  }
}

axis(1, xlab="Effect of mental load priming", cex.axis=1.5)
axis(2, at=c(6,5,4,3,2,1), labels=c('Political\ninterest', 'Political\nambition', 'Vote\nintention', 'Public\nparticipation', 'Private\nparticipation', 'Workplace\nadvancement'), las=2, cex.axis=1.5)

title(xlab="Effect of mental load priming", cex.lab=1.5)

# Close the tiff device
dev.off()


### Interpret results 

mean(dat$sinterest, na.rm=TRUE)
sd(dat$sinterest, na.rm=TRUE)

mod1a<-lm(slocal~treated, data=dat)
summary(mod1a)
mod1b<-lm(snational~treated, data=dat)
summary(mod1b)
mod1c<-lm(sinternational~treated, data=dat)
summary(mod1c)

mean(dat$work_scale1_2, na.rm=TRUE)
sd(dat$work_scale1_2, na.rm=TRUE)

####################################
### Figure 4
####################################

mod7<-lm(work_hrs~treated, data=dat)
summary(mod7)
mod8<-lm(work_hrs_partner~treated, data=dat)
summary(mod8)

results <- list(mod7, mod8)
fx <- sapply(results, function(x){unlist(summary(x)$coefficients["treated", 1])})
l95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 1])})
u95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 2])})
ys <- c(1,2)

tiff("Fig4.tiff", height = 6, width = 10, units = "in", res = 600, compression = "lzw")

par(mar=c(5.1,10,2.1,2.1))

plot(fx, ys, xlim=c(-3.5,1),ylim=c(0,3), ylab='', yaxt='n', type='n', xaxt='n', main="", xlab="")
abline(v=0)
for (i in 1:2){
  if(i %% 2 != 0){
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  }
  else {
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  }
}
axis(1, xlab="Effect of mental load priming", cex.axis=1.5)
axis(2,at=c(1,2), labels=c('Work\nhours','Partner\nwork\nhours'), las=2, cex.axis=1.5)
title(xlab="Effect of mental load priming", cex.lab=1.5)

dev.off()

### descriptive patterns, do men and women differ in preferred working hours of partner ? 

t.test(work_hrs_partner ~ woman, data=subset(dat, treated==0))
t.test(work_hrs_partner ~ woman, data=subset(dat, treated==1))

check<-lm(work_hrs_partner~woman, data=subset(dat, treated==0))
summary(check)

## how does it compare to preferences for own working hours? 

t.test(work_hrs ~ woman, data=subset(dat, treated==0))
t.test(work_hrs ~ woman, data=subset(dat, treated==1))

############################################
##### Figure 5
############################################

mod1a<-lm(sinterest~treated, data=subset(dat, woman==1))
summary(mod1a)
mod1b<-lm(sinterest~treated, data=subset(dat, woman==0))
summary(mod1b)
mod2a<-lm(spol_ambition_scale~treated, data=subset(dat, woman==1))
summary(mod2a)
mod2b<-lm(spol_ambition_scale~treated, data=subset(dat, woman==0))
summary(mod2b)
mod3a<-lm(spol_vote_scale~treated, data=subset(dat, woman==1))
summary(mod3a)
mod3b<-lm(spol_vote_scale~treated, data=subset(dat, woman==0))
summary(mod3b)
mod4a<-lm(spol_public~treated, data=subset(dat, woman==1))
summary(mod4a)
mod4b<-lm(spol_public~treated, data=subset(dat, woman==0))
summary(mod4b)
mod5a<-lm(spol_private~treated, data=subset(dat, woman==1))
summary(mod5a)
mod5b<-lm(spol_private~treated, data=subset(dat, woman==0))
summary(mod5b)
mod6a<-lm(work_scale1_2~treated, data=subset(dat, woman==1))
summary(mod6a)
mod6b<-lm(work_scale1_2~treated, data=subset(dat, woman==0))
summary(mod6b)

results <- list(mod1a,mod1b, mod2a, mod2b, mod3a, mod3b, mod4a, mod4b, mod5a, mod5b, mod6a, mod6b)

fx <- sapply(results, function(x){unlist(summary(x)$coefficients["treated", 1])})
l95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 1])})
u95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 2])})
ys <- c(6,5.75,5,4.75,4,3.75, 3,2.75, 2,1.75, 1, .75)

tiff("Fig5.tiff", height = 8, width = 10, units = "in", res = 600, compression = "lzw")
par(mar=c(7,10,2.1,2.1))

plot(fx, ys, xlim=c(-.15,.1), ylab='', yaxt='n', type='n', xaxt='n', main="", xlab="")
abline(v=0)
for (i in 1:12){
  if(i %% 2 != 0){
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  }
  else {
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1, lty='dashed')
    points(fx[i], ys[i], pch=17, col='black', cex=2)
  }
}
axis(1, xlab="Effect of mental load priming", cex.axis=1.5)
axis(2,at=c(6,5,4,3,2,1), labels=c('Political\ninterest','Political\nambition', 'Vote\nintention', 'Public\nparticipation', 'Private\nparticipation', 'Workplace\nadvancement'), las=2, cex.axis=1.5)
legend(x="bottom", inset = c(0, -0.2), ncol=1, bty='n', legend=c('Women', 'Men'), lty=c('solid','dashed'), lwd = 2, pch=c(16,17), cex=1.5,
text.width=c(0.0,0.05), x.intersp=0.5, horiz=TRUE, xpd=T)
title(xlab="Effect of mental load priming by gender", cex.lab=1.5) 

dev.off()

########################################
#### Figure 6
########################################

mod7a<-lm(work_hrs~treated, data=subset(dat, woman==1))
summary(mod7a)
mod7b<-lm(work_hrs~treated, data=subset(dat, woman==0))
summary(mod7b)
mod8a<-lm(work_hrs_partner~treated, data=subset(dat, woman==1))
summary(mod8a)
mod8b<-lm(work_hrs_partner~treated, data=subset(dat, woman==0))
summary(mod8b)

results <- list(mod7a, mod7b, mod8a, mod8b)
fx <- sapply(results, function(x){unlist(summary(x)$coefficients["treated", 1])})
l95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 1])})
u95 <- sapply(results, function(x){unlist(confint(x, level=0.95)["treated", 2])})
ys <- c(2, 1.75, 1, 0.75)

tiff("Fig6.tiff", height = 8, width = 10, units = "in", res = 600, compression = "lzw")
par(mar=c(7,10,2.1,2.1))

plot(fx, ys, xlim=c(-5,2),ylim=c(0,3), ylab='', yaxt='n', type='n', xaxt='n', main="", xlab="")
abline(v=0)
for (i in 1:4){
  if(i %% 2 != 0){
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1)
    points(fx[i], ys[i], pch=16, col='black', cex=2)
  }
  else {
    segments(l95[i],ys[i], u95[i],ys[i], col='black', lwd=1, lty='dashed')
    points(fx[i], ys[i], pch=17, col='black', cex=2)
  }
}
axis(1, xlab="Effect of mental load priming", cex.axis=1.5)
axis(2,at=c(2, 1), labels=c('Work\nhours','Partner\nwork\nhours'), las=2, cex.axis=1.5)
legend(x="bottom", inset = c(0, -0.23), ncol=1, bty='n', legend=c('Women', 'Men'), lty=c('solid','dashed'), lwd = 2, pch=c(16,17), cex=1.5,
text.width=c(0.0,1), x.intersp=2, horiz=TRUE, xpd=T)
title(xlab="Effect of mental load priming by gender", cex.lab=1.5) 

dev.off()

###################
## Open-ended text analysis 
###################

dat$open2_norespond<-0
dat$open2_norespond[dat$open2_count==0]<-1
table(dat$open2_norespond)
211/(791+211) ## 21% responded

dat$open2_respond<-1
dat$open2_respond[dat$open2_norespond==1]<-0
table(dat$open2_respond)
(dat$open2_respond)

res <- t.test(open2_respond ~ woman, data=dat)
res ## no significant difference by gender in likelihood of responding

## believe mental load negatively impacts politics / workplace advancement interest
## for this i will only use those responses that express a clear answer to the question we asked about this 

table(dat$open_2_coding)

dat$open_yes<-NA
dat$open_yes[dat$open_2_coding=="no"]<-0
dat$open_yes[dat$open_2_coding=="yes"]<-1
table(dat$open_yes)

res <- t.test(open_yes ~ woman, data=dat)
res ## yes, there is a significant gender difference in the responses. Although the majority of men and women who
## answered this question believe that it does. 81% women vs 57% men. 


###################
###################
## Appendix #######
###################
###################

############# Table A1

### Balance check table for Appendix (Table A1)
t1 <- t.test(woman~ treated, data=dat)
t2 <- t.test(highered ~ treated, data=dat)
t3 <- t.test(white~ treated, data=dat)
t4 <- t.test(Age~ treated, data=dat)
t5 <- t.test(full_time ~ treated, data=dat)
t6 <- t.test(age_youngest_1 ~ treated, data=dat)
t7 <- t.test(no_child2 ~ treated, data=dat)
t8 <- t.test(ml_share~ treated, data=dat)

tab <- map_df(list(t1, t2, t3, t4, t5, t6, t7, t8), tidy)

balance<-tab[c("estimate1", "estimate2", "estimate", "p.value")]
xtable(balance, type = "latex")

############## Table A2 

stargazer(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8,
          digits=2, dep.var.labels.include = T, no.space=T,
          omit.stat = c("f", "rsq", "ser"))


############### Table A3 mothers only

stargazer(mod1a, mod2a, mod3a, mod4a, mod5a, mod6a, mod7a, mod8a,
          digits=2, dep.var.labels.include = T, no.space=T,
          omit.stat = c("f", "rsq", "ser"))

############## Table A4 fathers only  

stargazer(mod1b, mod2b, mod3b, mod4b, mod5b, mod6b, mod7b, mod8b,
          digits=2, dep.var.labels.include = T, no.space=T,
          omit.stat = c("f", "rsq", "ser"))


############### Table A5 Controlling for negative emotions

mod1c<-lm(sinterest~treated + sneg_emotions, data=dat)
mod2c<-lm(spol_ambition_scale~treated + sneg_emotions, data=dat)
mod3c<-lm(spol_vote_scale~treated + sneg_emotions, data=dat)
mod4c<-lm(spol_public~treated + sneg_emotions, data=dat)
mod5c<-lm(spol_private~treated + sneg_emotions, data=dat)
mod6c<-lm(work_scale1_2~treated + sneg_emotions, data=dat)
mod7c<-lm(work_hrs~treated + sneg_emotions, data=dat)
mod8c<-lm(work_hrs_partner~treated + sneg_emotions, data=dat)

stargazer(mod1c, mod2c, mod3c, mod4c, mod5c, mod6c, mod7c, mod8c,
          digits=2, dep.var.labels.include = T, no.space=T,
          omit.stat = c("f", "rsq", "ser"))

### Among women only (not shown in text) 
#mod1<-lm(sinterest~treated + sneg_emotions, data=subset(dat, woman==1))
#mod2<-lm(spol_ambition_scale~treated + sneg_emotions, data=subset(dat, woman==1))
#mod3<-lm(spol_vote_scale~treated + sneg_emotions, data=subset(dat, woman==1))
#mod4<-lm(spol_public~treated + sneg_emotions, data=subset(dat, woman==1))
#mod5<-lm(spol_private~treated + sneg_emotions, data=subset(dat, woman==1))
#mod6<-lm(work_scale1_2~treated + sneg_emotions, data=subset(dat, woman==1))
#mod7<-lm(work_hrs~treated + sneg_emotions, data=subset(dat, woman==1))
#mod8<-lm(work_hrs_partner~treated + sneg_emotions, data=subset(dat, woman==1))


############# Table A6: Interaction with negative emotions 

mod1d<-lm(sinterest~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod2d<-lm(spol_ambition_scale~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod3d<-lm(spol_vote_scale~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod4d<-lm(spol_public~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod5d<-lm(spol_private~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod6d<-lm(work_scale1_2~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod7d<-lm(work_hrs~treated + sneg_emotions + treated:sneg_emotions, data=dat)
mod8d<-lm(work_hrs_partner~treated + sneg_emotions + treated:sneg_emotions, data=dat)


stargazer(mod1d, mod2d, mod3d, mod4d, mod5d, mod6d, mod7d, mod8d,
          digits=2, dep.var.labels.include = T, no.space=T,
          omit.stat = c("f", "rsq", "ser"))


### Among women (not shown in main text)
#mod1<-lm(sinterest~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod2<-lm(spol_ambition_scale~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod3<-lm(spol_vote_scale~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod4<-lm(spol_public~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod5<-lm(spol_private~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod6<-lm(work_scale1_2~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod7<-lm(work_hrs~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))
#mod8<-lm(work_hrs_partner~treated + sneg_emotions + treated:sneg_emotions, data=subset(dat, woman==1))


###########################################
###########################################
### Figure 2 = topic models ###############
###########################################
###########################################

library(ggplot2)
library(tidyr)
library(plyr)
library(dplyr)
library(haven)
library(stargazer)
library(xtable)
library(plm)
library(grid)
library(coefplot)
library(lubridate)
library(ggpubr)
library(modelsummary)
library(stringi)
library(quanteda)
library(tidyverse)
#install.packages("topicmodels")
library(topicmodels)
#install.packages("ldatuning")
library(ldatuning)
library(stm)
library(wordcloud)
library(reshape2)

theme_set(theme_bw())
set.seed(12345)

load("prolific_data_bjps_replication.RData")
head(dat)
is.factor(dat$woman)
dat$Gender<-NA
dat$Gender[dat$woman==0]<-"Fathers"
dat$Gender[dat$woman==1]<-"Mothers"
dat$Gender<-as.factor(dat$Gender)

vars<-c("Gender", "mental_load_open", "ml_share")
tp <-dat[,vars]
tp2<-na.omit(tp)

tokens <- tp2$mental_load_open %>%
  tokens(what = "word",
         remove_punct = TRUE,
         remove_numbers = TRUE,
         remove_url = TRUE) %>%
  tokens_tolower() %>%
  tokens_remove(stopwords("english"))

dfm <- dfm_trim(dfm(tokens), min_docfreq = 0.005, max_docfreq = 0.99, 
                docfreq_type = "prop", verbose = TRUE)

dfm_stm <- convert(dfm, to = "stm")

model <- stm(documents = dfm_stm$documents,
             vocab = dfm_stm$vocab, 
             K = 5,
             verbose = TRUE)

summary(model)

###############################
### Figure 2 in text 

tiff("Fig2.tiff", height = 6, width = 14, units = "in", res = 900, compression = "lzw")

plot(model, n=9,  topic.names = c("Cleaning: ", "Scheduling: ", "Financial / Home maintenance: ", "Anticipating needs: ", 
"Child care / school: "), xlim = c(0, .7))

dev.off()
###############################

### try to understand what each topic is about (used for labels)

tp3<-tp2[1:994, ]

thoughts1 <- findThoughts(model,texts=tp3$mental_load_open,
topics=1, n=3)$docs[[1]]
plotQuote(thoughts1) ##cleaning 

thoughts2 <- findThoughts(model,texts=tp3$mental_load_open,
topics=2, n=3)$docs[[1]]
plotQuote(thoughts2) ##scheduling 

thoughts3 <- findThoughts(model,texts=tp3$mental_load_open,
topics=3, n=3)$docs[[1]]
plotQuote(thoughts3) ##financial planning

thoughts4 <- findThoughts(model,texts=tp3$mental_load_open,
topics=4, n=3)$docs[[1]]
plotQuote(thoughts4)  ##anticipating household needs

thoughts5 <- findThoughts(model,texts=tp3$mental_load_open,
topics=5, n=3)$docs[[1]]
plotQuote(thoughts5) ## child care / school


model2 <- stm(documents = dfm_stm$documents,
             vocab = dfm_stm$vocab, 
             K = 5,
prevalence =~ Gender, content =~ Gender,
             max.em.its = 75, 
data = tp3,
init.type = "Spectral")

summary(model2)

thoughts1 <- findThoughts(model2,texts=tp3$mental_load_open,
topics=1, n=3)$docs[[1]]
plotQuote(thoughts1) ##cleaning 

thoughts2 <- findThoughts(model2,texts=tp3$mental_load_open,
topics=2, n=3)$docs[[1]]
plotQuote(thoughts2) ##scheduling 

thoughts3 <- findThoughts(model2,texts=tp3$mental_load_open,
topics=3, n=3)$docs[[1]]
plotQuote(thoughts3) ##financial /home maintenance

thoughts4 <- findThoughts(model2,texts=tp3$mental_load_open,
topics=4, n=3)$docs[[1]]
plotQuote(thoughts4)  ##anticipating household needs

thoughts5 <- findThoughts(model2,texts=tp3$mental_load_open,
topics=5, n=3)$docs[[1]]
plotQuote(thoughts5) ## child care / school


plot(model2, type="perspectives",  topics = c(1))
plot(model2, type="perspectives",  topics = c(2))
plot(model2, type="perspectives",  topics = c(3))
plot(model2, type="perspectives",  topics = c(4))
plot(model2, type="perspectives",  topics = c(5))

##### Gender differences 

prep <- estimateEffect(1:5 ~ Gender, model2,
 meta = tp3, uncertainty = "Global")

summary(prep, topics=1)
summary(prep, topics=2) ##mothers more likely 
summary(prep, topics=3) ##fathers more likely 
summary(prep, topics=4)
summary(prep, topics=5)





